/**
 * Computer Graphics Lecture WS 2009/2010 Ulm University
 * Creator: Manuel Finckh
 * Email:   manuel.finckh@uni-ulm.de
 */

#include "render.h"

float ray_prob = 2.0;

float n1 = 1.0;
float n2 = 1.33;

float actual_n = n1;

Render::Render(Scene *scene) : scene(scene) {
    accel = new BVH(scene->triangles, scene->num_tris);

    cam = scene->cam;
    ResX = cam->ResX;
    ResY = cam->ResY;

    sResX = SNAPX;
    sResY = SNAPY;

    accum_buffer = new int[ResX * ResY];

    nr_samples = 0;

    image = new Vec3[ResX * ResY];
    for (int i = 0; i < ResX * ResY; i++) image[i] = Vec3(0.0f, 0.0f, 0.0f);

    for (int i = 0; i < ResX * ResY; i++) accum_buffer[i] = 0;

    snapshot = new Vec3[sResX * sResY];
    for (int i = 0; i < sResX * sResY; i++) snapshot[i] = Vec3(0.0f, 0.0f, 0.0f);

    mtrand = new MTRand(1337);
    accum_index = 0;
}

Render::~Render() {
    delete accel;
    delete [] image;
    delete mtrand;
}

inline void Render::clearImage() {
    for (int i = 0; i < ResX * ResY; i++) {
        image[i] = Vec3(0.0f, 0.0f, 0.0f);
        accum_buffer[i] = 0;
    }
    nr_samples = 0;
}

float* Render::render(int shader, int snap) {
    float inv_accum, shrink;
    shrink_accum(inv_accum, shrink);

    actual = image;
    nResX = ResX;
    nResY = ResY;
    if (cam->moved) {
        for (int i = 0; i < nResX * nResY; i++)
            actual[i] = Vec3(0.0f, 0.0f, 0.0f);
        clearImage();
    }
    depth = 0;

    for (int y = 0; y < nResY; y++) {
        for (int x = 0; x < nResX; x++) {
            if (shader == 4) {
                int id;
                Vec3 light_point = scene->sampleRandomLight(id, mtrand->randExc(), mtrand->randExc(), mtrand->randExc());
                Vec3 dir2(0.0);
                //cout<<id<<endl;
                dir2 = sample_hemisphere(scene->triangles[id].getNormal());
                dir2.normalize();
                Ray ray(light_point, dir2, RAY_EPS, RAY_MAX);
                HitRec rec = accel->intersect(ray);
                // * scene->num_lights * scene->triangles[id].area()
                light_trace(ray, rec, scene->material[scene->mat_index[id]].color_e * scene->num_lights * scene->triangles[id].area());
                nr_samples++;

            } else {
                actual[x + y * nResX] *= shrink;
                float xc = x + mtrand->randExc(), yc = y + mtrand->randExc();

                Ray ray = cam->getRay(xc, yc);
                HitRec rec = accel->intersect(ray);
                Vec3 hp = cam->eye + rec.dist * ray.dir;
                //float x1, y1;
                //cam->getPoint(hp-cam->eye,x1,y1);
                //cout<<"x: ("<<x<<" , "<<x1<<")   y:("<<y<<" , "<<y1<<")"<<endl;
                Vec3 color(0.0f, 0.0f, 0.0f);
                //cout<<"x: "<<x<<"   y: "<<y<<endl;
                if (rec.id != -1) {
                    switch (shader) {
                        case 0: color = shade_debug(ray, rec);
                            break;
                        case 1: color = shade_simple(ray, rec);
                            break;
                        case 2: color = shade_path(ray, rec);
                            break;
                        case 3: color = shade_pathDL(ray, rec);
                            break;
                        case 5: color = shade_path_SSS(ray, rec);
                            break;
                        default: color = shade_simple(ray, rec);
                            break;
                    }
                } else if (scene->environment != NULL) {
                    color = scene->getEnvironment(ray.dir);
                    //ray.dir.printVec3();
                }
                actual[x + y * nResX] += color * inv_accum;
                actual_n = n1;
            }
        }
    }



    return (float*) actual;
}

Vec3 Render::shade_debug(Ray &ray, HitRec &rec) {
    //Vec3 normal = scene->triangles[rec.id].getNormal();
    Vec3 normal = scene->getShadingNormal(ray, rec.id);
    normal.abs();
    if (scene->material[scene->mat_index[rec.id]].tex != NULL) {
        // TODO
    }
    return normal;
}

Vec3 Render::shade_simple(Ray &ray, HitRec &rec) {
    //Vec3 normal = scene->triangles[rec.id].getNormal();
    Vec3 normal = scene->getShadingNormal(ray, rec.id);
    float cos = fabsf(normal * ray.dir);
    Vec3 color = cos * scene->material[scene->mat_index[rec.id]].color_d;
    if (scene->material[scene->mat_index[rec.id]].tex != NULL) {
        float alpha, beta, gamma;
        scene->triangles[rec.id].getBarycentric(ray, alpha, beta);
        gamma = 1.0f - alpha - beta;
        Vec2 coords = gamma * scene->uv[rec.id * 3] + alpha * scene->uv[rec.id * 3 + 1] + beta * scene->uv[rec.id * 3 + 2];

        color = Vec3::product(color, scene->material[scene->mat_index[rec.id]].tex->getColor(coords));
    }
    return color;
}

Vec3 Render::shade_mirror(Ray &ray, HitRec &rec) {
    Vec3 normal = scene->getShadingNormal(ray, rec.id);
    Vec3 mirror = scene->triangles[rec.id].reflect(ray, normal);
    Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
    Ray reflect_ray(hitpoint, mirror, RAY_EPS, RAY_MAX);
    HitRec r2 = accel->intersect(reflect_ray);

    Vec3 color(0.0f, 0.0f, 0.0f);
    if (r2.id != -1) {
        return shade_mirror(reflect_ray, r2);
    } else if (scene->environment != NULL) {
        color = scene->getEnvironment(reflect_ray.dir);
    }
    return color;

}

Vec3 Render::shade_path(Ray &ray, HitRec &rec) {

    Vec3 color_e(0.0f);
    rec.setDepth(rec.depth + 1);
    if (rec.id != -1 && rec.depth <= 5) {




        /*if(scene->material[scene->mat_index[rec.id]].name.compare("Glass") == 0){
                //cout<<"test"<<endl;
                scene->material[scene->mat_index[rec.id]].glass(dir2, ray.dir, normal,mtrand->randExc());
        }
        else*/
        //if(mtrand->randExc() < ray_prob){
        //Vec3 normal = scene->triangles[rec.id].getNormal();
        Vec3 normal = scene->getShadingNormal(ray, rec.id);
        if (scene->material[scene->mat_index[rec.id]].color_e.length_sqr() > 0.0001) {
            if (ray.dir * normal < 0.0)
                color_e = scene->material[scene->mat_index[rec.id]].color_e;
            else
                normal = -normal;

        }
        //Vec3 normal = scene->getShadingNormal(ray, rec.id);

        //normal.normalize();
        Vec3 color = scene->material[scene->mat_index[rec.id]].color_d;
        //Vec3 mirror = scene->triangles[rec.id].reflect(ray, normal);
        Vec3 hitpoint = ray.origin + rec.dist * ray.dir;



        Vec3 dir2(0.0);
        dir2 = sample_hemisphere(normal);

        dir2.normalize();
        Ray reflect_ray(hitpoint, dir2, RAY_EPS, RAY_MAX);
        HitRec r2 = accel->intersect(reflect_ray);
        r2.setDepth(rec.depth);
        // texture ?
        /*if ( scene->material[scene->mat_index[rec.id]].tex != NULL ) {
                float alpha, beta, gamma;
                scene->triangles[rec.id].getBarycentric(ray, alpha, beta); gamma = 1.0f - alpha - beta;
                Vec2 coords = gamma * scene->uv[rec.id * 3] + alpha * scene->uv[rec.id * 3 + 1] + beta * scene->uv[rec.id * 3 + 2];
	    
                color = Vec3::product(color, scene->material[scene->mat_index[rec.id]].tex->getColor(coords));
        }*/

        // fabs(reflect_ray.dir*normal)*
        return (color_e + Vec3::product(color, shade_path(reflect_ray, r2))); // / ray_prob;

        //}
        /*else if ( scene->environment != NULL ) {
                //color = scene->getEnvironment(ray.dir);
                color = Vec3(0.0);
        }	*/
    }
    return Vec3(0.0);


}

inline Vec3 Render::direct_illumination(Ray &ray, HitRec &rec) {
    Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
    int id;
    Vec3 light_dir = (scene->sampleRandomLight(id, mtrand->randExc(), mtrand->randExc(), mtrand->randExc()) - hitpoint);
    float dist = light_dir.length_sqr();
    light_dir.normalize();
    Vec3 light_normal = scene->triangles[id].getNormal();
    if (light_dir * light_normal < 0.0) {
        Ray light_ray(hitpoint, light_dir, RAY_EPS, RAY_MAX);
        HitRec light_rec = accel->intersect(light_ray);

        // not in shadow
        if (light_rec.id == id) {
            Vec3 radiance(0.0);
            //Vec3 normal = scene->triangles[rec.id].getNormal();
            Vec3 normal = scene->getShadingNormal(ray, rec.id);
            float geo_term = fabs(normal * light_dir) * fabs(light_normal * light_dir);
            geo_term /= dist;
            //cout<<light_dir.length_sqr()<<endl;
            if (scene->material[scene->mat_index[rec.id]].name.compare("SSS") == 0) {
                radiance = scene->material[scene->mat_index[id]].color_e * geo_term *
                        WardBRDF(-ray.dir, light_dir, normal) * scene->num_lights * scene->triangles[id].area();
            } else {
                radiance = scene->material[scene->mat_index[id]].color_e * geo_term *
                        BRDF(-ray.dir, light_dir, normal) * scene->num_lights * scene->triangles[id].area();
            }
            //** scene->num_lights * scene->triangles[id].area() ;
            //cout<<geo_term<<endl;
            //radiance.printVec3();
            //scene->material[scene->mat_index[id]].color_e.printVec3();
            return radiance;
        }
    }

    return Vec3(0.0);


}
// point + triangle id

inline Vec3 Render::direct_illumination(Vec3 & p, int &p_id) {
    Vec3 hitpoint = p;
    int id;
    Vec3 light_dir = (scene->sampleRandomLight(id, mtrand->randExc(), mtrand->randExc(), mtrand->randExc()) - hitpoint);
    float dist = light_dir.length_sqr();
    light_dir.normalize();
    Vec3 light_normal = scene->triangles[id].getNormal();
    if (light_dir * light_normal < 0.0) {
        Ray light_ray(hitpoint, light_dir, RAY_EPS, RAY_MAX);
        HitRec light_rec = accel->intersect(light_ray);

        // not in shadow
        if (light_rec.id == id) {
            Vec3 radiance(0.0);
            //Vec3 normal = scene->triangles[p_id].getNormal();
            Vec3 normal = scene->getShadingNormal(light_ray, p_id);
            float geo_term = fabs(normal * light_dir) * fabs(light_normal * light_dir);
            geo_term /= dist;
            //cout<<light_dir.length_sqr()<<endl;
            radiance = scene->material[scene->mat_index[id]].color_e * geo_term * (1 / M_PI) * scene->num_lights * scene->triangles[id].area();
            return radiance;
        }
    }

    return Vec3(0.0);


}

inline Vec3 Render::shade_pathDL(Ray &ray, HitRec &rec) {
    //cout<<rec.depth<<endl;
    Vec3 color(0.0f);
    rec.setDepth(rec.depth + 1);
    if (rec.id != -1 && rec.depth <= 10) {
        if (scene->material[scene->mat_index[rec.id]].color_e.length_sqr() > 0.0001) {
            if (rec.depth <= 1 && (scene->triangles[rec.id].getNormal() * ray.dir) < 0.0)
                return scene->material[scene->mat_index[rec.id]].color_e;
            else
                return Vec3(0.0);
        }
        //Vec3 normal = scene->triangles[rec.id].getNormal();
        Vec3 normal = scene->getShadingNormal(ray, rec.id);
        //color = scene->material[scene->mat_index[rec.id]].color_d;
        //color = Vec3(1.0);

        Vec3 hitpoint = ray.origin + rec.dist * ray.dir;

        Vec3 direct = direct_illumination(ray, rec);

        //return Vec3::product(color,fabs(ray.dir*normal));
        if (mtrand->randExc() < ray_prob) {
            Vec3 dir2(0.0);
            color = scene->material[scene->mat_index[rec.id]].color_d;

            if (scene->material[scene->mat_index[rec.id]].name.compare("Glass") == 0) {
                //cout<<"glas"<<endl;
                float R, T;
                //cout<<"ray * normal: "<<ray.dir * normal<<endl;
                if (ray.dir * normal > 0.0) {
                    fresnel(R, T, ray.dir, -normal, n1, n2);
                    //cout<<"von innen"<<endl;

                } else {
                    fresnel(R, T, ray.dir, normal, n2, n1);
                    //cout<<"von außen"<<endl;
                }
                //cout<<"R: "<<R<<"   T: "<<T<<endl;
                float r = mtrand->rand();
                //cout<<"r: "<<r<<endl;
                if (r < R)
                    dir2 = scene->triangles[rec.id].reflect(ray);

                else {
                    //cout<<"ref: ";
                    if (ray.dir * normal > 0.0)
                        refract(dir2, ray.dir, -normal, n1 / n2);
                    else
                        refract(dir2, ray.dir, normal, n2 / n1);

                    //cout<<dir2*normal<<endl;;
                }
                direct = Vec3(0.0);
                color = Vec3(1.0);

            } else {
                //cout<<"kein glas"<<endl;
                dir2 = sample_hemisphere(normal);
            }
            //dir2 = sample_hemisphere(normal);
            dir2.normalize();
            Ray reflect_ray(hitpoint, dir2, RAY_EPS * 10, RAY_MAX);
            HitRec r2 = accel->intersect(reflect_ray);
            r2.setDepth(rec.depth);


            //color.printVec3();
            //direct.printVec3();
            return Vec3::product(color, (direct + shade_pathDL(reflect_ray, r2)) / ray_prob);
        }


    } else if (scene->environment != NULL)
        //color = scene->getEnvironment(ray.dir);
        color = Vec3(0.0);

    return color;
}

inline float Render::rad2deg(float rad) {

    return rad * 180.0 / M_PI;
}

inline void Render::light_trace(Ray &ray, HitRec &rec, Vec3 radiance) {

    Vec3 color_e(0.0f);
    if (rec.id != -1) {
        Vec3 normal = scene->triangles[rec.id].getNormal();
        if (scene->material[scene->mat_index[rec.id]].color_e.length_sqr() > 0.0001) {
            if (ray.dir * normal < 0.0)
                color_e = scene->material[scene->mat_index[rec.id]].color_e;
            else
                normal = -normal;
        }
        Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
        Vec3 cam_dir = hitpoint - cam->eye;
        Ray shadow_ray(cam->eye, cam_dir, RAY_EPS, RAY_MAX);

        HitRec shadow_rec = accel->intersect(shadow_ray);

        // not in shadow

        if (shadow_rec.id == rec.id) {
            float _x, _y, _weight;
            cam->getPoint(cam_dir, _x, _y, _weight);

            //cout<<_x<<"    "<<_y<<endl;
            //cout<<_weight<<endl;

            if (_x >= 0.0 && _x < cam->ResX && _y >= 0.0 && _y < cam->ResY) {
                if (cam->w * cam_dir > 0.0) {
                    Vec3 color = scene->material[scene->mat_index[rec.id]].color_d;
                    float dist = cam_dir.length_sqr();
                    cam_dir.normalize();
                    float geo_term = fabs(normal * cam_dir) * fabs(cam_dir * cam->w);
                    geo_term /= dist;

                    accumulate((int) _x, (int) _y, Vec3::product((color), radiance * geo_term * BRDF(-ray.dir, cam_dir, normal) * _weight));
                }

            }


        }
        //if(mtrand->randExc() < ray_prob){

        Vec3 dir2(0.0);
        Vec3 color = scene->material[scene->mat_index[rec.id]].color_d;
        dir2 = sample_hemisphere(normal);
        dir2.normalize();
        Ray reflect_ray(hitpoint, dir2, RAY_EPS, RAY_MAX);
        HitRec r2 = accel->intersect(reflect_ray);
        float cos_l = fabs(dir2 * normal);
        light_trace(reflect_ray, r2, (Vec3::product(color, radiance)));
        //}
    }

}

inline Vec3 Render::shade_path_iterative(Ray &ray, HitRec &rec) {
    Vec3 color = Vec3(0.0f);
    Vec3 path_weight(1.0f);
    int depth = 0;

    while (true) {
        Vec3 point = ray.origin + rec.dist * ray.dir;
        Vec3 s_normal = scene->getShadingNormal(ray, rec.id);
        Vec3 g_normal = scene->triangles[rec.id].getNormal();
        Material *mat = &scene->material[scene->mat_index[rec.id]];

        float cos_g = g_normal * ray.dir;
        if (cos_g > 0.0f) g_normal = -g_normal;
        if (s_normal * g_normal < 0.0f) s_normal = -s_normal;
        float cos_s = -(s_normal * ray.dir);
        if (cos_s < 0.0f) break;
        if (cos_g < 0.0f) {
            color += Vec3::product(mat->color_e, path_weight);
        }
        if (++depth > 5) break;

        path_weight = Vec3::product(path_weight, mat->color_d);

        Vec3 dir_new;
        Material::diffuse(dir_new, s_normal, mtrand->rand(), mtrand->rand());
        ray.origin = point;
        ray.dir = dir_new;
        ray.tmin = RAY_EPS;
        ray.tmax = FLT_MAX;

        rec = accel->intersect(ray);
        if (rec.id == -1) break;
    }
    if (scene->environment != NULL) {
        color += Vec3::product(path_weight, scene->getEnvironment(ray.dir));
    }
    return color;
}

inline float Render::rand_exponential(const float &lambda) {
    return (-1.0 / lambda) * log(mtrand->rand());
}

inline float Render::schlick(const Vec3 &in, const Vec3 &n) {


    float theta = fabs(in * n);
    float r0 = ((n1 - n2) / (n1 + n2)) * ((n1 - n2) / (n1 + n2));

    return (r0 + (1.0 - r0) * pow((1.0 - cos(theta)), 5.0));
}

float phase = 1.0 / (4 * M_PI);

float albedo = 0.5;

float sigma_s = 2.0;
float sigma_a = 0.00035;
float sigma_t = sigma_s + sigma_a;
float sigma_tr = sqrtf(3.0 * sigma_a * sigma_t);

inline float Render::dipole(Vec3 x0, Vec3 normal, Vec3 &x1, int &id) {

    float Fdr = 1.440 / (n2 * n2) + 0.710 / n2 + 0.668 + 0.0636 * n2;

    float albedo_r = sigma_s / sigma_t;

    float zr = 1.0 / sigma_t;
    float A = (1.0 + Fdr) / (1.0 - Fdr);
    float D = 1.0 / (3.0 * sigma_t);
    float zv = zr + 4.0 * A * D;

    Vec3 u, v, w;
    Vec3::onb(u, v, w, normal);
    int z = 1;
    HitRec v1;
    Ray upper;
    float d;
    int cnt = 10;
    do {
        d = rand_exponential(sigma_tr);
        //cout<<d<<endl;
        Vec3 l = mtrand->rand() * u + mtrand->rand() * v;
        l.normalize();
        x1 = x0 + d * l;

        upper = Ray(x1 + z * normal, -normal, RAY_EPS, RAY_MAX);
        //Ray lower(x0 -  normal, x1 - (x0 -  normal), RAY_EPS, RAY_MAX);
        v1 = accel->intersect(upper);
    } while (v1.id == -1 && cnt-- > 0);

    id = v1.id;
    x1 = upper.origin + v1.dist * upper.dir;
    d = (x1 - x0).length();
    //cout<<d<<endl;

    float dr = sqrtf(zr * zr + d * d);
    float dv = sqrtf(zv * zv + d * d);
    float term1 = (sigma_tr * dr + 1) * exp(-sigma_tr * dr) / (sigma_t * pow(dr, 3.0));
    float term2 = zv * (sigma_tr * dv + 1) * exp(-sigma_tr * dv) / (sigma_t * pow(dv, 3.0));

    return (albedo_r / (4.0 * M_PI) * (term1 + term2));



}

inline Vec3 Render::shade_path_SSS(Ray &ray, HitRec &rec) {
    //cout<<rec.depth<<endl;
    Vec3 color(0.0f);
    rec.setDepth(rec.depth + 1);
    if (rec.id != -1 && rec.depth <= 5) {
        if (scene->material[scene->mat_index[rec.id]].color_e.length_sqr() > 0.0001) {
            if (rec.depth <= 1 && (scene->triangles[rec.id].getNormal() * ray.dir) < 0.0)
                return scene->material[scene->mat_index[rec.id]].color_e;
            else
                return Vec3(0.0);
        }
        //Vec3 normal = scene->triangles[rec.id].getNormal();
        Vec3 normal = scene->getShadingNormal(ray, rec.id);


        Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
        Vec3 SSS(0.0);


        if (scene->material[scene->mat_index[rec.id]].name.compare("SSS") == 0) {
            int z = 5;

            float s0 = rand_exponential(sigma_t);
            Vec3 refr(0.0);

            HitRec inner;
            Ray bounce_ray;
            refract(bounce_ray.dir, ray.dir, normal, n1 / n2);
            inner.dist = s0;
            bounce_ray.origin = hitpoint;
            while (z--) {
                Vec3 P = bounce_ray.origin + inner.dist * bounce_ray.dir;
                int id;
                Vec3 light_P = scene->sampleRandomLight(id, mtrand->randExc(), mtrand->randExc(), mtrand->randExc());

                Ray light_ray(light_P, P - light_P, RAY_EPS, RAY_MAX);
                HitRec light_rec = accel->intersect(light_ray);


                if (light_rec.id >= scene->material[scene->mat_index[rec.id]].start && light_rec.id <= scene->material[scene->mat_index[rec.id]].end) {
                    //normal = scene->triangles[light_rec.id].getNormal();
                    normal = scene->getShadingNormal(ray, light_rec.id);
                    Vec3 in_P = light_ray.origin + light_rec.dist * light_ray.dir;
                    float s1 = (in_P - P).length();
                    // nachher
                    float c_i = fabs(light_ray.dir * normal);
                    float s1_prime = s1 * (c_i / (sqrtf(1.0 - (1.0 / n2)*(1.0 / n2) * (1.0 - (c_i * c_i)))));

                    //Vec3 direct(0.0);

                    float t1, t2, r;
                    fresnel(r, t1, ray.dir, normal, n2, n1);
                    fresnel(r, t2, light_ray.dir, normal, n2, n1);

                    //float geo_term = 1.0; //fabs(normal * ray.dir) / fabs(normal * light_ray.dir);
                    Vec3 light_normal = scene->triangles[id].getNormal();
                    float geo_term = fabs(light_normal * light_ray.dir);
                    geo_term /= (light_rec.dist * light_rec.dist);
                    float sigma_tc = sigma_t; // + (fabs(normal * ray.dir) / fabs(normal * light_ray.dir)) * sigma_t;
                    //cout<<light_dir.length_sqr()<<endl;

                    SSS = scene->material[scene->mat_index[id]].color_d * scene->num_lights * geo_term * scene->triangles[id].area();
                    SSS *= (sigma_s * t1 * t2 * phase / sigma_tc) * exp(-s1_prime * sigma_t) * exp(-s0 * sigma_t) / M_PI;
                    color += SSS;
                    s0 = rand_exponential(sigma_t);


                    bounce_ray = Ray(P, sample_sphere(), 0.0, RAY_MAX);
                    inner = accel->intersect(bounce_ray);
                    if (inner.id != -1) {
                        if (inner.dist < s0) {
                            z = 0;
                            // später vllt: refract
                            hitpoint = bounce_ray.origin + inner.dist * bounce_ray.dir;
                            normal = scene->triangles[inner.id].getNormal();
                            rec = inner;
                            ray = bounce_ray;
                        } else {
                            inner.dist = s0;
                        }

                    }
                }
            }
        }

        Vec3 dir2(0.0);
        dir2 = sample_hemisphere(normal);
        dir2.normalize();
        // get direct illumination
        Vec3 direct = direct_illumination(ray, rec);

        // dipol zeug
        if (scene->material[scene->mat_index[rec.id]].name.compare("SSS") == 0) {
            float t1, t2, r;
            //fresnel(r,t1,ray.dir,normal,n1,n2);
            //fresnel(r,t2,dir2,normal,n1,n2);
            t1 = t2 = 1.0;
            int x_id;
            Vec3 x1;
            float Rd = dipole(hitpoint, normal, x1, x_id);
            Vec3 x_direct = direct_illumination(hitpoint, x_id);
            direct += Rd * x_direct * t1 * t2 * (1.0 / M_PI) * (1.0 / sigma_tr);
        }/**/


        //if(mtrand->randExc() < ray_prob){

        Ray reflect_ray(hitpoint, dir2, RAY_EPS, RAY_MAX);
        HitRec r2 = accel->intersect(reflect_ray);
        r2.setDepth(rec.depth);
        color = scene->material[scene->mat_index[rec.id]].color_d;
        return Vec3::product(color, (direct + shade_path_SSS(reflect_ray, r2)));
        //}


    } else if (scene->environment != NULL)
        color = Vec3(0.0);

    return color;
}


// n1 = actual medium, n2 = adjacent one

inline void Render::fresnel(float &R, float &T, const Vec3 &in, const Vec3 &n, float n1, float n2) {

    float n0 = n1 / n2;
    float c_i = fabs(in * n);
    float s_t = n0 * n0 * (1.0 - c_i * c_i);
    float c_t = sqrtf(1.0 - s_t);
    //cout<<"critical: "<<rad2deg(asin(n0))<<"  winkel: "<<rad2deg(acos(c_i))<<endl;
    if (s_t > 1.0) {
        R = 1.0;
        T = 0.0;
        return;
    }

    float R_r = (n1 * c_i - n2 * c_t) / (n1 * c_i + n2 * c_t);
    float R_t = (n2 * c_i - n1 * c_t) / (n2 * c_i + n1 * c_t);

    R = 0.5 * (R_r * R_r + R_t * R_t);
    T = 1.0 - R;



}

inline void Render::refract(Vec3 &out, const Vec3 &in, const Vec3 &n, float n0) {


    float c1 = n * (-in);
    float c2 = sqrt(1 - (n0 * n0) * (1 - (c1 * c1)));
    //if(c1 > 0)
    out = n0 * in + (n0 * c1 - c2) * n;
    //else
    //	out = n0 * in + (n0 * c1 + c2) * n;

}

inline void Render::shrink_accum(float &inv_accum, float &shrink) {
    if (!cam->moved) {
        accum_index++;
        inv_accum = 1.0f / (float) accum_index;
        shrink = (accum_index - 1) * inv_accum;
    } else {
        accum_index = 1;
        inv_accum = 1.0f;
        shrink = 0.0f;
        for (int i = 0; i < ResX * ResY; i++)
            accum_buffer[i] = 1;
        clearImage();

    }
}

inline void Render::accumulate(const int &x, const int &y, const Vec3 &color) {

    /*
    Vec3 t = actual[x + y * ResX];
    int N = accum_buffer[x + y * ResX];
    actual[x + y * ResX] = (N * t + color) / (N+1);
    accum_buffer[x + y * ResX]++;
     */

    Vec3 t = actual[x + y * ResX];
    int N = nr_samples;
    actual[x + y * ResX] = (N * t + color) / (N + 1);


}

inline Vec3 Render::sample_hemisphere(Vec3 normal) {
    /*float theta = mtrand->randExc(2.0 * M_PI);
    float phi = mtrand->randExc(M_PI / 2.0);
    float r = cos(phi);
    Vec3 tmp(r*cos(theta),r*sin(theta),sin(phi));
    tmp.normalize();*/
    float u = mtrand->randExc();
    float v = mtrand->randExc();
    Vec3 u1, v1, w1, out;
    Vec3::onb(u1, v1, w1, normal);

    const float phi = 2.0f * M_PI * u;
    const float sqrtu = sqrtf(v);
    const float x = sqrtu * cosf(phi);
    const float y = sqrtu * sinf(phi);
    const float z = sqrtf(1.0f - v);

    out = x * u1 + y * v1 + z * w1;
    return out;
}

inline Vec3 Render::sample_sphere() {

    float u = mtrand->randExc();
    float v = mtrand->randExc();
    Vec3 out(0.0);


    const float phi = 2.0f * M_PI * u;
    const float theta = M_PI * v;
    out.x = sin(theta) * cos(phi);
    out.y = sin(theta) * sin(phi);
    out.z = cos(theta);
    return out;
}


float alpha = 0.071, rho_s = 0.07, rho_d = 0.60;
float alpha_s = alpha * alpha;
// Ward BRDF

inline float Render::WardBRDF(Vec3 in, Vec3 out, Vec3 normal) {

    in.normalize();
    out.normalize();
    Vec3 half = in + out;
    half.normalize();
    Vec3 x1, y1, z1;
    Vec3::onb(x1, y1, z1, normal);
    float gamma = tan(acos(half * normal));
    float s = rho_s / (4.0 * alpha_s * M_PI * sqrtf(fabs((in * normal) * (out * normal))));
    //std::cout<<"BRDF(s): "<<4.0 * alpha * alpha *  (in * normal) * (out * normal)<<std::endl;
    //float exponent = ((half * x1) * (half * x1) / alpha_s + (half * y1) * (half * y1)/ alpha_s)  / ((half * normal)*(half * normal));
    float exponent = gamma * gamma / alpha_s;
    //std::cout<<"BRDF: "<<((rho_d + s * exp(-exponent)))<<std::endl;
    return (rho_d / M_PI + s * exp(-exponent));
    //return 1.0/M_PI;
    //return 1.0;

}

inline float Render::BRDF(Vec3 in, Vec3 out, Vec3 normal) {

    return 1.0 / M_PI;

}





